/*
 * Copyright (C) 2010 - Alexandru Gagniuc - <mr.nuke.me@gmail.com>
 * This file is part of ElectroMag.
 *
 * ElectroMag is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * ElectroMag is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 *  along with ElectroMag.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "Electrostatics.h"
#include "Electrodynamics.h"
#include "Magnetics.h"
#include "Data Structures.h"
#include "X-Compat/HPC Timing.h"


/**
 * \brief Confines a particle to a rectangular box in the limits boxMin, boxMax
 */
template<class T>
inline void BoundToBox (
    Vector3<T>& position,
    Vector3<T>& velocity,
    const Vector3<T> boxMin,
    const Vector3<T> boxMax
)
{
    if ( position.x > boxMax.x )
    {
        position.x = 2 * boxMax.x - position.x;
        velocity.x = -velocity.x;
    }
    if ( position.y > boxMax.y )
    {
        position.y = 2 * boxMax.y - position.y;
        velocity.y = -velocity.y;
    }
    if ( position.z > boxMax.z )
    {
        position.z = 2 * boxMax.z - position.z;
        velocity.z = -velocity.z;
    }

    if ( position.x < boxMin.x )
    {
        position.x = 2 * boxMin.x - position.x;
        velocity.x = -velocity.x;
    }
    if ( position.y < boxMin.y )
    {
        position.y = 2 * boxMin.y - position.y;
        velocity.y = -velocity.y;
    }
    if ( position.z < boxMin.z )
    {
        position.z = 2 * boxMin.z - position.z;
        velocity.z = -velocity.z;
    }
}

/**
 * \brief Calculates the electro-magnetic and gravitational fields generated by
 * \brief a dynamic particle, and adds the result the existing field
 */
template<class T>
inline void emgFieldsAdd (
    const electro::dynamicPointCharge<T>& srcpart,
    const Vector3<T>& position,
    Vector3<T>& electroField,
    Vector3<T>& magneticField,
    Vector3<T>& gravitField
)
{
    Vector3<T> rInvSq = vec3InverseSquare (
        vec3 ( position,srcpart.staticProp.position ) );
    //vec3InverseSquareFLOP (10 FLOP)

    // Electric field contribution
    electroField += electro::PartFieldOp (
        srcpart.staticProp.magnitude, rInvSq );
    // (5 FLOP)

    // Magnetic field contribution
    magneticField += magnetic::PartFieldOp (
        srcpart.velocity, srcpart.staticProp.magnitude, rInvSq );
    // (13 FLOP)

    // TODO: gravitational field contribution
    // gravitField += gravitFieldOp();
    // (5 FLOP)

    // Total: 33 FLOP / jstep (28 w/o grav);
}

/**
 * \brief executes one Verlet integration step on an
 * \brief electro-magneto-gravitational particle system
 */
template<class T>
int emgVerletStep (
    Array<electro::dynamicPointCharge<T> >& dynCharges,
    Array<electro::pointCharge<T> >& pointCharges,
    Array<Vector3<T> >& verletAccel,
    Vector3<T> minBound,
    Vector3<T> maxBound,
    const T timeStep,
    const size_t n,
    perfPacket& perfData
)
{
    //TODO: Quick parameter check;
    if ( !n )
        return 1;

    //get the size of the computation
    size_t p = pointCharges.GetSize();

    const T dt = timeStep;
    const T dt2 = dt / 2;

    // Big loop
    for ( size_t i = 0; i < n; i++ ) //one particle at a time
    {
        electro::dynamicPointCharge<T> ipart = dynCharges[i];
        Vector3<T> ipartAcc = verletAccel[i];

        //TODO: v(t + dt/2) = v(t) + a(t) * dt/2;
        ipart.velocity += ipartAcc * dt2;
        // 6 FLOP

        //TODO: r(t + dt) = r(t) + v(t + dt/2) * dt;
        ipart.staticProp.position += ipart.velocity * dt;
        // 6 FLOP

        //TODO: Force computation F( r(t + dt) );
        Vector3<T> electroField = {0, 0, 0};
        for ( size_t jstat = 0; jstat < p; jstat++ )
        {
            electroField += electro::PartField ( pointCharges[jstat],
                                                 ipart.staticProp.position );
        }
        Vector3<T> magneticField = {0, 0, 0};
        Vector3<T> gravitField = {0, 0, 0};
        // Use to loops to avoid calculating particle interacting with itself
        for ( size_t j = 0; j < i; j++ )
        {
            emgFieldsAdd ( dynCharges[j], ipart.staticProp.position,
                           electroField, magneticField, gravitField );
            // Total: 33 FLOP / jstep (28 w/o grav);
        }
        for ( size_t j = i + 1; j < n; j++ )
        {
            emgFieldsAdd ( dynCharges[j], ipart.staticProp.position,
                           electroField, magneticField, gravitField );
            // Total: 33 FLOP / jstep (28 w/o grav);
        }

        ipartAcc = ipart.staticProp.magnitude * electroField +      // 3 FLOP
                   magnetic::Force ( ipart, magneticField );        // 10 FLOP
        // gravitForce;                                     // 3 FLOP
        // 16 FLOP (13 w/o grav);

        //TODO: a(t + dt) = F(cur)/m(cur);
        ipartAcc /= ipart.mass;
        // 3 FLOP

        //TODO: v(t + dt) = v(t + dt/t) + a(t + dt) * dt/2;
        ipart.velocity += ipartAcc * dt2;
        // 6 FLOP

        verletAccel[i] = ipartAcc;

        BoundToBox ( ipart.staticProp.position,
                     ipart.velocity,
                     minBound, maxBound );
    }
    // 37 FLOP + (n - 1) * 33 FLOP) / istep

    size_t totalFLOP = ( 37 + ( n - 1 ) * 33 ) * n;
    return 0;
}

